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Abstract. - We develop a theory for the full counting statistics (FCS) for a class of nano- 
electromechanical systems (NEMS), describable by a Markovian generalized master equation. 
The theory is applied to two specific examples of current interest: vibrating Ceo molecules 
and quantum shuttles. We report a numerical evaluation of the first three cumulants for the 
C6o-setup; for the quantum shuttle we use the third cumulant to substantiate that the giant en- 
hancement in noise observed at the shuttling transition is due to a slow switching between two 
competing conduction channels. Especially the last example illustrates the power of the FCS. 



Introduction. - The full counting statistics (FCS) of charge transport in mesoscopic sys- 
tems is an active topic of recent research [1-5]. Calculation and measurement of the whole 
probability distribution of transmitted charge is motivated by the fact that FCS provides more 
information about a particular system than just the mean current or current noise which are 
the first two cumulants of the large-time asymptotics of the probability distribution. Very 
recently, a measurement of the third cumulant, which quantifies the skewness of the distri- 
bution, was reported [6]. The detailed nature of charge transport in nanoelectromechanical 
systems (NEMS), another modern field in mesoscopics, poses many challenges both to exper- 
iments and theory, and the computation of FCS for NEMS is an important task that needs 
to be addressed. The first steps were taken recently with a calculation of FCS for a driven, 
classical shuttle [7] . 

In this Letter, we present a theory for the evaluation of cumulants in a wide class of 
NEMS encompassing the majority of systems considered thus far, namely those which can be 
described by a Markovian generalized master equation (GME). The current cumulants turn 
out to be fully determined by an extremal eigenvalue of the system evolution superopera- 
tor (Liouvillcan) in analogy with previous studies [4,5]. Their evaluation is, however, more 
complicated since in NEMS there are generally many relevant states which need to be taken 
into account. We solve the problem by formulating a systematic perturbation theory, and 
using this derive explicit formulas for the first three cumulants. The method is illustrated by 
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two examples of NEMS - the Ceo-experiment [8] and the quantum shuttle [9-11]. To test the 
method we calculate the first three cumulants for the model of the Cgo-setup from [12]. In case 
the oscillator is equilibrated the cumulants can be calculated alternatively using P(S)-theory 
which gives the same results. For the quantum shuttle we use the third cumulant to substan- 
tiate that the giant enhancement of the current noise in the transition region [10] is caused 
by a slowly fluctuating amplitude of the shuttle resulting in a slow switching between two 
current channels, i.e. tunneling and shuttling. This part complements [7], which considered 
a fixed driving amplitude, and [13] describing a related phenomenon in a different model. 

Theory. - We consider a nanoelectromechanical system with discrete energy levels elec- 
tronically coupled to two leads and mechanically coupled to a generic heat bath providing 
dissipation. The system is described by the reduced density operator p(t), which we assume 
evolves according to the Markovian GME 

m = m)- a) 

The Liouvillean C describes the dynamics of the system and we assume that the system tends 
exponentially to a stationary state p stat . This implies that the Liouvillean, which is a non- 
hermitian operator, has a single eigenvalue equal to zero with p stat being the corresponding 
(normalized and unique) right eigenvector which we denote by |0)) [16]. The corresponding left 
eigenvector is the identity operator 1, denoted by ((0|, and we have ((0|0)) = Tr(F/3 stat ) = 1. 
The pair of eigenvectors allows us to define the projectors V = |0))((0| and Q = 1 — V obeying 
the relations VC = CP = and QCQ = C. We also introduce the pseudoinverse of the 
Liouvillean 1Z = QC~ 1 Q 1 which is well-defined, since the inversion is performed only in the 
subspace spanned by Q, where £ is regular. The assumption of exponential decay to the 
stationary state is equivalent to the spectrum of C in the subspace spanned by Q having a 
finite negative real part. 

In order to evaluate the FCS of the system, i. e. the probability P n (t) of n electrons being 
collected in, say, the right lead in the time span t, we resolve the density operator p(t) and the 
GME with respect to n. The GME is a continuity equation for the probability (charge) and, 
therefore, we can identify terms corresponding to charge transfer processes between the system 
and the right lead. Specifically, we introduce the superoperator 2 + of the particle current of 
electrons tunneling from the system to the right lead, and the corresponding superoperator 
T~ of the reverse process, where electrons tunnel from the right lead to the system. In terms 
of these superoperators the n-resolved GME can be written as 

3 (n) (*) = {C-1+ -l-)p {n) {t)+l+p {n - 1 \t)+l-p {n+1 \t) (2) 

with n = . . . , —1,0, 1, From the n-resolved density operator we can obtain, at least in 

principle, the complete probability distribution P n (t) = Tr [p(")(t)]. 

It is practical first to evaluate the cumulant generating function S(t, x) defined as 

oo 

e s(t,x) = P n {t)e mx . (3) 

n— — oo 

From S(t,x) we then find the m'th cumulant of the charge distribution (we take e = 1) 
by taking the m'th derivative with respect to the counting field x at x = 0, i.e. ({n m ))(t) = 
dfi X ) m 1*=°' an< ^ lrom tne knowledge of all cumulants we can reconstruct P n (t). The cumulants 
of the current in the stationary limit t — > oo are given by the time derivative of the charge 
cumulants, i.e. ((i™)) = ^(( nm ))(i)| ■ The first two current cumulants give the average 
current and the zero-frequency current noise, respectively. 
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Using p(")(t) we may express S(t,x) as e s{t ^ = Trf^^L^ P (n) {t)e inx ] = Tr[F(t,x)], 
where we have introduced the auxiliary operator F(t, \) whose equation of motion follows 
from the n-resolved GME, 

^F(t, X ) = [C-(1- e^)l+ (1 - e-^)I-}F(t, X ) = C x F(t 7 X ) (4) 

with the formal solution F(t, x) — e c * t F{§ 1 x). We assume adiabatic evolution of the spectrum 
of L x with increasing \j *- e - there is a unique eigenvalue A™ m of C x associated with the 
projector V x which develops from the zero eigenvalue of C and which is the closest to zero for 
small enough \- The rest of the spectrum still has finite negative real part which ensures the 
damping of its contribution for large times. Thus, we have 

e S(t, x ) = ((o| e ***|F(0, x )» - e A x in *((0|^ x |F(0,x))) = e^ t+c ^ for t - oo, (5) 

where C™* depends on the initial state of the system. However, the current cumulants in 
the stationary state do not depend on the initial conditions, but are totally determined by 
A™ m in full analogy with previous studies [4,5]. For NEMS in general C x is of very large 
dimensions, and the numerical evaluation of higher order derivatives of A™ ln may become a 
formidable numerical problem. In order to circumvent this problem we determine A™ m using 
Rayleigh-Schrodinger perturbation theory for C x = C + £' , treating C' x as the perturbation. 
Since the Liouvillean is not hermitian, we cannot assume that it has a spectral decomposition 
in terms of its eigenvectors, and one cannot use directly the standard formulas. However, it 
is possible to formulate the perturbation theory exclusively in terms of the projectors V, Q 
and the pseudoinverse 1Z. As in standard Rayleigh-Schrodinger perturbation theory the first 
order correction is given by the average of the perturbation with respect to the unperturbed 
eigenstatc. Taking the derivative of the first order correction with respect to ix and letting 
X — > we find 

((i)) = «6| X |0», (6) 

where I = T + — 2~ . As expected the first cumulant equals the average current. For the 
second cumulant, i.e. the zero- frequency current noise, one finds 

(U 2 )) = ((6|J|0))-2((0|^X|0)), (7) 

where J = T + +T~. In the high bias limit (where ((0|X _ |0)} = 0, since backward tunneling 
is blocked) this expression yields the result previously derived in [16]. The expression for the 
third cumulant 

{{I 3 )) = ((0|X |0)) - 3({0\111J + JK1 10)) - 6((0\JTZ(nir - IK)I |0)) (8) 

is the main result of this section, and below we evaluate it for two specific cases. Higher order 
cumulants can be obtained in the same manner by calculating the corresponding higher order 
corrections. 

Model 1: The Cgo experiment. - An experiment with a NEMS that has received much 
attention is the measurement of the TV-curves of a vibrating Ceo-molecule [8] . The experiment 
has been modelled in several papers [12, 14, 15] using a model which will also be employed 
here. Calculations of TV-curves [12, 14] have been found to be in good agreement with the 
experiment, and the current noise has been predicted [15]. We calculate the third cumulant 
for this setup by applying our method to the model as described in [12]. 
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Fig. 1 - Results for the C6o-setup. First three cumulants as function of the bias V. The parameters, 
which correspond to fig. 3 of [12], are E c = lOmeV, Kuj = 5meV, hV {0) = 1 ^eV, k B T = 0.15 meV, 
K = O.lwo, and ci = (dashed), 0.8 (dotted), 1.5 (full), c 2 = 0.005. 



In the model of [12] both the coupling to the leads and to the heat bath are treated in 
the weak coupling approximation which reduces the full GME to an ordinary Pauli master 
equation for the probabilities of occupation of the individual eigenstates of the system: 



dP„ 



dt 



[Wi 



P, 



m' .a' ,1 



(9) 



where m,a,l denote the (extra) charge on the molecule (m = 0,1), the spin, and the vi- 
brational state, respectively, while s indicates whether an electron tunnelled from/to the left 
(s = — 1) or right (s = 1) lead. P m ,a,i is the probability of being in the eigenstate labelled by 
the subindices. Bath-mediated transitions between different vibrational states are given by 
the thermal rates 



Wi +M = Wi^ l+1 e- h ^' kBT = K- 



l + l 



e hu> /k B T _ i •■ 

where ujq is the natural oscillator frequency. The charge transfer rates are 
r (s) „fm,,,„ ^t.si,,, ,o „,„ seV 



(10) 



l,o-,Z'<-0,0,Z 



T^\(l'\e^-^\l)\ 2 f(E c + ^- + nw {l' - I - 7 2 )), 



seV 



0,0,1^1,0,1* = r(°)|(^( at - a >|0| 2 [l - f(E c + — + huo(l' I - 7 2 ))], 



(11) 



where / is the Fermi function, is the bare tunneling rate, Ec is the charging energy 
difference, and V is the symmetrically applied bias. The quantity 7 describes the bias- 
dependence of the electric field at the position of the molecule, and is assumed to have the 
form 7 = ci + t|^C2 [12]. Here we do not consider the case where the rates depend on the 
position of the molecule, although this can easily be included. 

The current superoperators are identified from the expression for the stationary current 



T-stat \ ^ "p(t) pstat 

1 ~ |_ 0,0,1' <~l,cr,l r l,<J,l\ 



,1.1' 



E 

a,l,V 



-p(l) pstat 
[ L 1, a, I' ^0,0,1^0,0,, 



(12) 



«0|X+|0» 



«0|X-|0» 
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Fig. 2 - The quantum shuttle consists of a nanosized grain moving in a harmonic potential between 
two leads. A high bias between the leads drives electrons through the grain. 

Here |0)) is a diagonal density matrix containing the stationary probabilities P^ 3 ^ ; ■ Due to this 
diagonal form of the density matrices the relevant superoperators needed for the cumulants 
can be represented by matrices of dimension 2N x 2N (N is the number of vibrational modes) 
which makes the calculation of the cumulants numerically straightforward. 

In flg.EJwe show the bias-dependence of the first three cumulants for parameters corre- 
sponding to fig. 3 of [12]. Since <C K, the oscillator is in equilibrium and, therefore, the 
FCS of the model can also be calculated from a simple two-level model with 4 rates given by 
P(_E)-theory. We have verified that the semi-analytical results coincide with numerics (not 
shown) which we view as a non-trivial test of our method. In non-equilibrium cases there are 
no simple alternatives to the numerics. We demonstrate the full power of the method in the 
second example. 

Model 2: The Quantum Shuttle. We consider the model of a quantum shuttle used 
in [9-11]. The system consists of an oscillating nanoscopic grain coupled to two leads (fig. 
EJ. In the strong Coulomb blockade regime the grain effectively has just one electronic level. 
The oscillations of the grain are treated fully quantum mechanically, and damping of the 
oscillations is due to a surrounding heat bath. As in [10] we consider the n-resolved system 
density matrices p\™\t), i = 0, 1, where n is the number of electrons that have tunneled into 

the right lead in the time span t. In the high bias limit n is nonnegative, and the p^(t) 
evolve according to the n-resolved GME 

&>(t) ^[^osc - eK,pW(t)]+£ dafflp ^(t) - ^L{e%,p^(t)} + r L e-ip$(t)e-i, 

(13) 

with n = 0, 1, . . . and p\ x X \t) = 0. Here the commutators describe the coherent evolution of 
the charged (pn) or empty (poo) shuttle which is modelled by a quantum mechanical harmonic 
oscillator of mass m and frequency ui. The electric field between the leads is denoted E. The 
terms proportional to T^m describe transfer processes from the left/to the right lead with 
hopping amplitudes that depend exponentially on the position ?, where A is the electron 
tunneling length. The mechanical damping of the oscillator is described by the damping 
kernel (here T = 0) C damp p = -%j-[x, {p, p}] - ^[x, [x,p\] [9, 10]. We identify the current 
superoperators from (fT^I : X + p — Tne* |0)(l|/5|l)(0|e* , I~ = 0. 

In [9, 10] it was found that the quantum shuttle exhibits a crossover from tunneling to 
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shuttling when the damping, starting above a certain threshold value, is decreased. This 
transition is clearly recorded both in the current [9] and the zero-frequency current noise [10]. 
The FCS in the tunneling and shuttling limit is to a first approximation captured by the re- 
sults for the zero amplitude (with appropriately renormalized rates) and the large (shuttling) 
amplitude of a driven shuttle [7], respectively. When approaching the semi-classical regime, a 
giant enhancement of the noise was found in the transition region. This behavior was tenta- 
tively attributed to amplitude fluctuations in the spirit of [13], however, a more quantitative 
explanation has been missing. The phase space representation of the stationary state of the 
shuttle in the transition region indicated that shuttling and tunneling processes coexist [10] 
leading to the conjecture that the giant noise enhancement is caused by switching between 
two current channels^) (tunneling and shuttling) induced by infrequent jumps between two 
discrete values of the shuttle amplitude. Very recently the FCS of such bistable systems has 
been studied [18], and it was found that the first three cumulants are (assuming that the 
individual channels are noiseless) 

((/» = ^t + ItF^s , (14a) 

t T^S + 1 S^T 

«/ 2 )} = 2(/ s - /t) 2 ,/ 5 ^ 1 ^ 5 , 3 , (14b) 

(Is^t + It^s)' 3 

((/ 3 )) = 6(i s - jT) 3 r ^Trr^ s (r T s - r s ^ T ) _ 



(rs^T + Tt 



-s 



Here Jg/T denote the current associated with the shuttling/tunneling channcl( 2 ), while Tt^s 
is the transition rate from the shuttling to the tunneling channel and Tg^T is the rate of the 
reverse transition. 

The rates can be evaluated analytically in the limit A — > oo and E — > [19]. However, 
for other parameters an alternative approach is needed. First, one evaluates the first three 
cumulants numerically following the theory presented above. Then, one calculates the rates 
from the first two cumulants (eqs. I|14all4b|l ) and compares the numerically calculated third 
cumulant from eq. JHJ) with the one obtained from eq. I|14c(l . The numerical calculation of the 
cumulants is only possible by a truncation of the oscillator Hilbert space. By retaining the N 
lowest oscillator states the (non-sparse) matrix representations of the superoperators entering 
eqs. (KilSt are of dimension 2N 2 x 2N 2 , which for the required values of N ~ 100 leaves us 
with non-trivial numerical matrix problems. However, using the iterative methods described 
in [16] the cumulants can be evaluated numerically. 

If fig. we show the 7-dependence of the first three cumulants for A = 1.5xo, where xq = 
yjh/muj. The first cumulant, the current, shows the transition from the tunneling to shuttling 
current with decreasing damping. The transition is also evident from the second cumulant, 
the zero-frequency current noise, which exhibits a giant enhancement in the transition region, 
before dropping to very low values in the shuttling region. Together with the numerical results 
for the third cumulant we show the analytic expression (eq. p4c|l ) with rates extracted from 
the first two cumulants. As can be seen the two data sets coincide, which we take as evidence 
that the quantum shuttle in the transition region indeed behaves as a bistable system for 
which the FCS is known [18]. When approaching the deep quantum regime, A ~ xo (not 
shown), the transition from tunneling to shuttling is smeared out and the distinct current 
channels cease to exist. In this limit the bistable system model is not valid. 

( 1 ) Such a behavior, referred to as the 'whistle' effect, was first reported in [17]. 

( 2 ) 7s = <V2tt, It = , with f R = r R e R / m ^^ 2 e 2 ^/" 1 " 2 *, f L = Y L e n l m ^ [10,19]. 
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Fig. 3 - Results for the quantum shuttle. First three cumulants as function of the damping 7. The 
parameters are Tl = Fr = O.OIoj, A = l.5xo, d = eE/mui 2 = 0.5so, xq — wh/mu. The full lines 
indicate numerical results, while the circles indicate the analytic expression for the third cumulant 
assuming that the shuttle in the transition region effectively behaves as a bistable system. 

Conclusion. We have presented a method for computation of the FCS for typical 
nanoelectromechanical systems and applied it to two specific models. For the Ceo-setup with 
equilibrated oscillator wc have calculated the first three cumulants and explained the results 
in terms of a simple two-level model. For the quantum shuttle we have used the first three 
cumulants as evidence that the shuttle in the transition region behaves as a bistable system. 
This example clearly illustrates the usefulness of the FCS in probing a microscopic system. 
Here we have only shown explicit expressions for the first three cumulants, our method, 
however, can be extended to the calculation of cumulants of any order, and we believe that 
the method has a broad range of applicability. 

* * * 
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manuscript and fruitful discussions. 
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